module probability

contains

    subroutine probx(signal, size, min, max, numbins, distx)
        implicit none

        integer, intent(in) :: size, numbins
        real, intent(in), dimension(0:size) :: signal
        real, intent(in) :: min, max
    
        real, intent(out), dimension(0:numbins) :: distx

        integer :: i, actsize
        real :: shift, binsize
    

        if(min<0) then
            shift = abs(min)
        else
            shift = -min
        endif

        actsize = 0

        binsize = (max-min)/numbins
        
        do i = 0, size
            if(signal(i) > min .and. signal(i) < max) then    
                distx( int(floor( (shift+signal(i)) / binsize)) ) = distx( int(floor( (shift+signal(i)) / binsize)) ) + 1.0
                actsize = actsize + 1  
            end if
        enddo

        do i = 0, numbins
            distx(i) = distx(i)/real(actsize)
            if(distx(i) < 0.0) then
                distx(i) = 0.0
            end if
        enddo
        

    end subroutine ! probx


    subroutine probxy(signala, mina, maxa, numbinsa, signalb, minb, maxb, numbinsb, distxy, size)
        implicit none

        integer, intent(in) :: size, numbinsa, numbinsb
        real, intent(in), dimension(0:size) :: signala, signalb
        real, intent(in) :: mina, maxa, minb, maxb
    
        real, intent(out), dimension(0:numbinsa, 0:numbinsb) :: distxy

        integer :: i, j, actsize
        real :: shifta, shiftb, binsizea, binsizeb

        if(mina<0) then
            shifta = abs(mina)
        else
            shifta = -mina
        endif
        
        if(minb<0) then
            shiftb = abs(minb)
        else
            shiftb = -minb
        endif

        actsize = 0

        binsizea = (maxa-mina)/numbinsa
        binsizeb = (maxb-minb)/numbinsb
        
        do i = 0, size
            if(signala(i) > mina .and. signala(i) < maxa .and. signalb(i) > minb .and. signalb(i) < maxb) then    
                distxy( int(floor( (shifta+signala(i)) / binsizea)), int(floor( (shiftb+signalb(i)) / binsizeb)) ) = &
distxy( int(floor( (shifta+signala(i)) / binsizea)), int(floor( (shiftb+signalb(i)) / binsizeb)) ) + 1.0

                actsize = actsize + 1  
            end if
        enddo

        do i = 0, numbinsa
            do j = 0, numbinsb
                distxy(i,j) = distxy(i,j)/real(actsize)
                if(distxy(i,j) < 0.0) then
                    distxy(i,j) = 0.0
                end if
            enddo
        enddo


    end subroutine ! probxy


    subroutine probxyz(signala, mina, maxa, numbinsa, signalb, minb, maxb, numbinsb, signalc, minc, maxc, numbinsc, distxyz, size)
        implicit none

        integer, intent(in) :: size, numbinsa, numbinsb, numbinsc
        real, intent(in), dimension(0:size) :: signala, signalb, signalc
        real, intent(in) :: mina, maxa, minb, maxb, minc, maxc
    
        real, intent(out), dimension(0:numbinsa, 0:numbinsb, 0:numbinsc) :: distxyz

        integer :: i, j, k, actsize
        real :: shifta, shiftb, shiftc, binsizea, binsizeb, binsizec

        
        if(mina<0) then
            shifta = abs(mina)
        else
            shifta = -mina
        endif
        
        if(minb<0) then
            shiftb = abs(minb)
        else
            shiftb = -minb
        endif
    
        if(minc<0) then
            shiftc = abs(minc)
        else
            shiftc = -minc
        endif

        actsize = 0

        binsizea = (maxa-mina)/numbinsa
        binsizeb = (maxb-minb)/numbinsb
        binsizec = (maxc-minc)/numbinsc
        
        do i = 0, size
            if(signala(i) > mina .and. signala(i) < maxa .and. signalb(i) > minb .and. &
signalb(i) < maxb .and. signalc(i) > minc .and. signalc(i) < maxc) then    
                distxyz( int(floor( (shifta+signala(i)) / binsizea)), int(floor( &
(shiftb+signalb(i)) / binsizeb)), int(floor( (shiftc+signalc(i)) / binsizec))) = &
distxyz( int(floor( (shifta+signala(i)) / binsizea)), int(floor( (shiftb+signalb(i)) / binsizeb)), &
int(floor( (shiftc+signalc(i)) / binsizec))) + 1.0
                actsize = actsize + 1  
            end if
        enddo

        do i = 0, numbinsa
            do j = 0, numbinsb
                do k = 0, numbinsc
                    distxyz(i,j,k) = distxyz(i,j,k)/real(actsize)
                    if(distxyz(i,j,k) < 0.0) then
                        distxyz(i,j,k) = 0.0
                    end if
                enddo
            enddo
        enddo

    end subroutine ! probxyz


end module
